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Abstract 

A variety of finite element schemes has been used in the numerical 
approximation of compressible flows particularly in underwater acoustics. In 
many instances instabilities have been generated due to the lack of mass 
conservation. In this paper we develop new two- and three-dimensional 

elements which avoid these problems. 
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INTRODUCTION 


It has been known for several years that many finite element formulations 
of fluid flows have potentially serious instabilities (see, e.g., [l]-[8]). A 
great deal of attention has been given to the incompressible case, and a 
variety of techniques for dealing with these instabilities has been 
developed. These include reduced integration [9], artificial diffusion [10], 
penalty methods [11], and the development of special stable elements [12]- 
[14], to cite only a few references in a very large literature. The special 
elements have proven to be of great practical value since they typically can 
be implemented by using standard elements in a Galerkin formulation, and then 
applying various filtering techniques to the pressures computed from the 
latter. 

The compressible case is similar in the sense that elements which are 
unstable for incompressible flows are also unstable in the compressible case. 
This has been analyzed mathematically and shown to be equivalent to the 
failure of the Babuska-Brezzi angle condition between the finite spaces used 
in the formulation [3]-[4] . On the other hand, the compressible case is 
dissimilar in that the angle condition is not sufficient for stability. Put 
differently, there are elements which perform quite satisfactorily for 
incompressible flows but which exhibit a number of instabilities in the 
compressible case. These are usually of the secondary nature; i.e., they do 
not destroy the entire calculation, but they do lead to slower rates of 
convergence than normal and, hence, less accurate and less efficient 
approximations . 

The key to these has been shown to be related to mass conservation [1], 
[4], [7]* In the compressible case it is necessary to have mass conservation 
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(i.e., a divergence free velocity field) only in an appropriate averaged 
sense. For example, if linear elements are used to represent the fluid 
velocity u^, typically one has 

) div u, = 0 (1.1) 
T 

for each subdivision T defining the grid. This is different from 

div = 0 (1.2) 

holding at each point, the latter being equivalent to exact mass conservation. 

In the steady compressible case, the relevant vector field is the mass 
flow u h , and the assertion is that without (1.2) serious errors will occur in 
the approximation that do not exist in the incompressible case. 

Fortunately, there are a number of mass conserving elements that have been 
developed in the two-dimensional case. The most widely used are the second- 
order accurate union jack or criss-cross element [3], and the first-order 
linear element developed by Thomas [14]. The goal of this paper is to develop 
analogous elements for three spatial dimensions. The analog of the Thomas 
element is given in Section 2, and the analog of the criss-cross element is 
given in Section 4. In Section 3 we give an apparently new element, which 
like the criss-cross element, is second-order accurate, but (like the Thomas 
element) does not have restrictive grid regularity conditions. 

To describe these elements we consider the specific case of steady 
potential flow. We let denote the potential and p the density. The 
governing equations in the flow region ft are: 
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div p grad $ = 

0 in 

Q y 

(1.3) 

♦ = 

<f> r on 

r D • 

(1.4) 

grad 4>*ri = 

v on 

n 

r N ‘ 

(1.5) 

On one part of the boundary of 

ft, namely 

Tp, the potential 

<|> r is given 


while on the other, T N> the normal velocity v n is given. The density p 
and velocity grad <j> are related through Bernoulli's equation. 

Galerkin formulations typically work in terms of the mass flow 

iu = p grad (f>, ( 1.6) 

since jump conditions across a shock are equivalent to the continuity of _u, 
and thus are an intrinsic part of the formulation. Entropy conditions ruling 
out expansion shocks, on the other hand, are introduced through modifications 
of the density p. The formulation consists of a finite element space in 

which the mass flow is represented and a finite space for the 

potential. The parameter h > 0 denotes an generic mesh spacing. One seeks 
a in Vj^ and a (j)^ in such that 

u, »v h . . 

J + } <J>. div v = J <(>_ v, *n (1.7) 

p 8 " 

J div u^ ty* 1 = 0. 


( 1 . 8 ) 
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The approximate mass flow 


is required to satisfy 


u *n = 
— n — 


u 

n 


on 


N’ 


(1.9) 


where u n is the given normal mass flow, and (1.7) holds for all v h in V h 
whose normal components are zero on 1^. The equation (1.8), on the other 
hand, holds for all in S^, the boundary condition (1.4) being natural in 

this formulation. Once a basis for and S h has been selected, (1.7)- 

(1.8) reduce to a system of nonlinear equations, the nonlinearity coming from 
the functional dependence of the density p on the mass flow u, . 

In the next three sections we shall display pairs V^, which satisfy 
exact mass conservations (1.2). In each of these cases this property is a 
direct consequence of the following inclusion property [3]: 


s h ■ div I v h' ' 


( 1 . 10 ) 


That is, each element i|> is a divergence div v|* of some v^ 1 in V^, and 
conversely. Indeed, suppose (1.10) were true for the pair S^, V^. Then if 
u, is the mass flow arising from (1.7)-(1.8) we have 


div 


% 


€ S. 


Thus letting i|; = div u^ in (1.8) we get 


} (div u, ) 2 = 0, 
0 


and hence (1.2) holds. 
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The inclusion property (1.10) is the primary vehicle in the paper for 
constructing mass conserving elements. The next two sections is taken as 

a suitable space of piecewise constant functions, and is constructed so 

that (1.10) holds. In the last section the converse is used. Here V h is 
taken as a suitable space of piecewise linear functions, and S h is defined 
by (1.10). 

Each of the spaces introduced in the next sections also satisfies the 
Babuska-Brezzi angle condition. The proofs will be omitted since they are 
very close to the two-dimensional proofs and quite technical in nature. 


2. A FIRST-ORDER ELEMENT 

We suppose for simplicity that 0 is a polynomial region in E? , and we 
subdivide it into tetrahedral elements, T 1 ,«»*,T n . The space S h consists 
of all piecewise constant functions on this grid. Thus has n degrees 
of freedom and a local basis standard manner, 


i.e.. 


( 1 if j = * 

( 0 if j = l . 


( 2 . 1 ) 


2 . 1 . 


To define V h we use the reference tetrahedron T ref shown in Figure 


- 6 - 



Figure 2 - l! T ref 

In particular, we consider incomplete linear functions in T re £ having the 
form 

v = a + 0 ? . (2.2) 


The four degrees of freedom a, 3 are uniquely determined by specifying the 
normal component v»v at the centroid m_ of each of the four faces of 
T ref* T * ie s P ace v h will consist of images of these functions. 

More precisely, let Tj be a tetrahedron in the grid, and consider the 
affine mapping 



( 2 . 3 ) 


Then functions in have the form in Tj : 
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= a + 35. 


(2.4) 


The parameters a, 3 are determined by specifying v«v at the centroid of 
each face. 

Suppose {Fp*",Fj is the set of all faces and let be the centroid 

of Fj with being a unit normal to Fj . Then has N degrees of 

freedom, and any in is uniquely determined by specifying 


vW 


( j = l » * * * >N) . 


(2.5) 


A basis £]. » * * * dlft for is obtained by requiring 


VVh * 5 ' 


KT 


if j - * 


if j * z 


( 2 . 6 ) 


where | F j | is the area of F j . Observe that since ^ is linear on Fj 


l V-t * 

b Z 


if j = Z 
if j * Z. 


(2.7) 


Moreover, <{>j is identically zero in any T^ which does not contain Fj as a 
face; i.e,, > * * * ».£ N is a locally defined basis. 

Functions ^ in are not necessarily continuous, but they are square 
integrable, and have a square integrable divergence. In fact, for each 
tetrahedron T 


dlv = tft l dlv ^ " W L *■ - 


in T, 


( 2 . 8 ) 
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where (T) denotes the volume of T and _v is the outer normal. This is 
the only regularity required on the variational principle (1.7)-(1.8). These 
elements, however, cannot be used for the incompressible case since continuity 
is required there. 

We assert that the inclusion property (1.10) holds for this pair of 
spaces. Since v^ is linear, div v^ is a piecewise constant and hence is 
in S h . We must show that every function ij> h has this form. To see this we 
first select the function v_ satisfying 

div v = t|/ h in 0. (2.9) 


Using the basis let 


N 




* i u 


j = l F 




( 2 . 10 ) 


Observe that for any T 


dlV ^'Tg-{ dlV 3i " TT7T It 


V 9T„ 


( 2 . 11 ) 


But using (2.7) and (2.10) we have 


9T 


} v^*_v =■ J _v*_v = J div v . 


9T 


( 2 . 12 ) 


In light of (2.9) we therefore have 


dlv ■ WJ J div ^ ■ TO- 1 *» 


h 


V T„ 


(2.13) 
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in T . Since \|», is an arbitrary function in S^, it follows that (1.10) 
x n 

must hold. 

Incidentally, the Babuska-Brezzi angle condition requires that, in 
addition to (2.13), Vjj must satisfy 

Hv^o + Udiv <_ CI^Iq (2.14) 

where 

lv« 0 - {} |v| 2 } 1/2 , (2.15) 

A 

and 0 < C < 00 is uniformly bounded independent of the grid. This can be 

verified directly from ( 2 . 10 ) once explicit formulas are obtained for the 

basis 6 , ,••• ,d>„. This calculation is similar to the one in 114], and it 
— 1 

does require that the grid be quasi- regular in the sense that the ratios of 
the maximum side length to the diameter of each T be uniformly bounded. 

Functions in are capable of only first-order accuracy. The reason 

for this can be seen from (2.2). The first component vj of v_ does not 
contain terms involving an< ^ £ 3 * Hence, v^ € is not a complete 

piecewise linear function. 

Since contains only piecewise constant functions, it too is capable 

of only first-order accuracy. Thus, if and ^ are the approximations 

obtained from (1.7)-(1.8), then 

llu - u^Iq < Chllgrad uH Q (2.16) 

(!((> - <f> h #Q < ChUgrad $Hq, (2.17) 


where h is the maximum diameter of tetrahedrons T^. 
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3. A SECOND-ORDER ELEMENT 

In many applications it is important to have second-order approximations 
to the mass flow and velocities. To achieve this we use the same grid as in 
Section 2, but with a larger space V h . The idea is to increase the size of 

n = {Vb 5 div ^ = 0}, (3.1) 

so that the inclusion property (1.10) is not affected. 

Indeed, we retain the same notation introduced in Section 2, and as 

before, let denote the space of piecewise constant functions. To define 

V h we let denote t ^ le faces of the tetrahedrons with 

being the associated centroids. With each face we select a normal v 

^ “J 

along with two independent tangential directions xf^ . A function v. 

— J — J — h 

in Vb> if and only if, 

(i) Vj, is a linear polynomial in each T j . 

(ii) Vb continuous at each centroid m^ . 

In the previous section, functions v^ were required only to have a 
continuous normal component at the centroids ra^ (j=l , • • • ,N) . Here all three 
components of vj^ are continuous. Moreover, the representation of Vb in 
each Tj is in terms of a complete linear polynomial. 

As in the last section, ^ in v h are not necessarily continuous 
everywhere in ft. There are, however, square integrable with square 

integrable divergences. Moreover, 

\ 

"''it.- ITT 


(3.2) 
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holds for each tetrahedron T. 

Observe that has 3N degrees of freedom, and a locally defined basis 

,j^*\i=0,l,2) can be selected so that for j«l,»««,N we 


ij 0) -VV - 



if j - i 
if j * 1 


(3.3) 


‘T^^On^) = 0 i = l * * * * »N, t=l,2, 


(3.4) 


and for i=l,2 


' 0 


all Z 


(3.5) 


.(i) (t), x 

J * * I 0 


if i = t and j = Z 


otherwise . 


(3.6) 


Note that since ^ is a H- near function on the faces we have 


J 4 0) Tk 

s z 


l 

0 


if j - i 
if j * Z 


(3.7) 


J = 0 all i,t (3.8) 

S * j 

and for i=l,2 



all l 


(3.9) 


1 


- 12 - 


S, *S“4» ■ { 


if j = l and i = t 
otherwise . 


Thus any _v^ in can be written 


v* - ? a ^M 0> ♦ ! I u 2h - i ' i, k' 1) . 

j=l Sj ^ J 2 j=l j=l s 2 2 


(3.10) 


(3.11) 


It follows from (3.2) and (3*9) that the second sura in (3.11) is divergence 
free* It therefore does not play a role in the inclusion property (1.10), but 
it does contain the extra terras (missing in the space defined in Section 2) 
giving second-order accuracy* 

To verify (1*10) for this case we select a function v satisfying (2*9), 
and let 

\ = I (J Z'lj k i (0> » (3.12) 

j=l S, J 2 


(i.e., the divergence free part is omitted). The equations (2.11 )— ( 2.13) 
remain valid, and hence as before we have 


div v^ = in T. 


(3.13) 


Assuming the grid is quasi-regular in the sense cited in Section 2, one 
can directly verify the angle condition introduced in [3]. A consequence is 
that second-order accuracy, i.e., 


"JL-Vo i 


< Ch II D ul, 


(3.14) 
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is obtained for the mass flow. Since S h consists only of piecewise constant 
elements, first-order accuracy (i.e., (2.17)) is the best one we can obtain to 
the potential <J>. However, the results in [3] indicate that <J>^ is a 

second-order accurate approximation to the average 

* 

in each tetrahedron T. 


4. THE THREE-DIMENSIONAL CRISS-CROSS ELEMENT 

Let 0 be a rectangular domain in IE?, divided into ra^ m 2 m 3 cubes 
where a = (i, j ,k )(0 _< i < , 0 <^j <m 2 » 0 ^ k < m^), by a uniform grid 

of mesh spacing h. We denote the corners of the cubes by centers 

by “i+ l / 2 J+ l / 2 ,k+ I /2 a " d the “ nteC ° f eaCh taCe by 3 l+ 1/2 , j+ 1/2 ,k' 

-i.j+‘/2 M l h or -1+ V 2 ,1 ,k+ V2 • 

Each cube is partitioned into 24 congruent tetrahedrons by considering 

divisions along the four diagonals of each cube, the two diagonals of each 

face of the cube and the three lines joining the centers of opposite faces of 

the cube [Figure 4.1] • We number the faces of each cube from 1 to 6 as shown 

in Figure 4.2. Each of these six faces of the cube will contain the bases of 

r s 

four tetrahedrons. We can now identify any tetrahedron in the grid by T * 

AW 

where a = (i,j,k) gives the cube in which the tetrahedron is contained, 

1 r 6 specifies the face on which the base is located, and 1 <, s < 4 

specifies the tetrahedron completely. 

Let be the space of all I? valued continuous functions over ft 

which are piecewise linear for each tetrahedron. We take our approximate 
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space to be that subspace of in which the nodes at the centers of 

the faces of each cube are condensed out. This is done by specifying the 

value at such a node to be the average of the values at the four vertices of 

3 

the face. More precisely, for v, = £ v^ e 0 ^^ e n. , 

h £=1 h 


1 1 


V “i+ V 2 ,j+ V 2 ,k } V i+ V 2 ,j+ V 2 ,k 1/4 q l Q V i+p,j+q,k ’ 


it =1,2,3. (4.1) 


The number v*(m i+ ^ ^ ^ ) and v*^ j+ ^ >k+ ^ ) will be similar 

averages • 

The dimensions of V h is equal to 3N, where N, the number of 

unconstrained nodes is given by (mj + l)(m 2 + lKmg + 1) + mj m 2 m 3 . Let n 

be the set of all nodes (constrained or unconstrained) . For every node m . 

—a’ 

A 

define a function <f> € V, by 

T a h 3 

W - S a 6 - 1 1£ s ■ £ 

/v A/ /V A/ 

=0 if a * 

h 3. 0 

Then any function v. n in V k can be represented as v n = \ ( £ v <j> )e - . 

£=1 a€ r\ ~ ~ 

For a function in V^, the coefficients for constrained nodes are replaced by 
suitable averages. 


( 1 ) 


±1 >® 2 ’-3 


for an orthonormal basis for I? 
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Let be the set of functions that are constant over each tetrahedron. 

<J» r,s will denote the value of d> in T r,S . Define Su to be the set of 
Y a Y a n 

rv rw 

/V 

<f> € such that for any a, 

(i) '•’a' 1 + ^a’ 3 “ ^a’ 2 + ^a’ 4 (r-l,*««,6) (4.2) 

/v /v cv 

4 

(ii) l + <f>^ ,S ) = C a> a constant (r = 1,2,3) (4.3) 

S=1 ~ ~ ~ 

With the choice of and V^> we have 33 div[V^] as proved in the 

appendix. 

Moreover, our proof also shows the Babuska-Brezzi condition is satisfied 
so that for every <(> € S^, we can find a w € satisfying 

div w = <f> in H 

nw« 0 < 

where 0 < C < 00 is bounded independently of the grid. 

The second-order accurate space of functions that are linear over each 

B is easily seen to be a subspace of V^. Similarly, the space of functions 

constant over each B is a subset of Sl, 

n n 
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APPEND IX 

We now prove the following theorem: 

Theefea: For the space and mentioned in Section 4, we have 

S h - div[V h ] . 

Moreover, for any <{> e S^, we can find a w € satisfying 


div w = j> in fl 


hwHq < c h < f> n — 1 . 


(A. 1) 


Proof : Let <j> = div v, where v £ V^. Referring to Figure A.l and 

1 s 

denoting A * by <|> s , we have 


Similarly, 


$1 - 


1 1 _ 2 2 2 „ 3 0 3 

V A - V 0 . 2V B ~ V A • V 0 . 2V 0 ■ 2v B 


, l l 1 2 2,3,3 

, 2s B - V 0 • V E ^ V E - V 0 , 2, c - 2v B 

<j, = r- + r + u 


1 1 2,2 „ 2 „3 „ 3 
t V D " V E . V D + V E " 2v B a 2v C ~ 2v B 
*3 h + h + K 


^A “ 


1 .. 1 0 1 2 2 0 3 0 3 

V D + V A " 2v B . V D ' V A . 2, C - 2v B 


From this, we see that 


♦l + *3 3 *2 + *4 = 


1^1 1 1 
V A + V D " V 0 " V E 


2 2 2 2 
V D + V E ' V 0 - V A 


+ 4 


, 3 3. 

(V C ' V 


(A. 2) 


This proves (4.2) for the function <|> in the case r = 1. The other cases 

can be verified similarly. We now prove (4.3) for <f>. We use <|> and <(>' 

s s 

1 S 

to denote <j> a * respectively. By calculations similar to the ones above, we 
obtain 


l 

8 


l 

S=1 


(♦s + +P 



!) ~ (vj + v^ 1 + v^ + v* 
4h 


0 


+ 




1 + V_ 1 + V„ 1 + V 


D 0 - ( 


3,3 
V A V 0 


+ V E + 


3 n 

V 


4h 


D v 1 D v^ + D, v^ = div, v 
x y 3 h 


(A. 3) 


3 3 

We have used averaging expressions for v B and v g in the above. By the 
symmetry of the expression for div^ v > it follows that the same result must 
be obtained for r = 2 and r = 3, i.e., (4.3) holds. Hence, <|> € S^, i.e., 
div[V h ] C S h . Now let <|> £ S h . We will find v h £ V h satisfying (A.3), 


which we rewrite as 
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Here, the left-hand side stands for the expression 


1 1 1 1 1 

J_ y y (fv 1 _ v 1 1 + fv 4 - v 4 ) 

L L U i+l,j+p,k+q i,j+p,k+q J ^ i+l,j+l,k+p i+q,j,k+p J 


4h 


q=0 p=0 


+ (v J - v 3 )\ 

L i+P.j+q.k+1 i+pj+q.k-’* 


D v 
x 


4.1 + D v 

i.j.k y 


4 4 1 + D v' 

i.j.k z 


i.j.k* 


(A. 5) 


The right-hand side equals 


4 J 

h J B 


div v. 


i.j .k 


Using (A. 2), we can obtain the value of at the central node as 


, 3 _ h r l,2 . .1,4 i 1 r If 1 1 'I 

i+ V 2 , j+ V 2 ,k+ V 2 A ^i» j »k ^i.j.k-l ” 4 ^ V i+l,j+p,k “ V i,j+p,k' 


( 2 _ 2 ^ / 3 ,3 

^ V i+P,j+l,k V i+p,j,k J ^ V i+p> j ,k V i+p,j+l,k^* 


(A. 6) 


Similar formulas hold for the other components of v_ at the central node of 

B . We shall use (A. 4) to determine v at the corner nodes and then define 
a 

v_ at central nodes by (A. 6). 
potential (« 1+ 1/ 2 l/ 2 , k+ 1/ 2 )• 

satisfying: 


To solve (A. 4), we introduce a discrete 
constant over each box ®i j k an d 
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1 1 


11*. 

i.j.k “ Ah q ^ 0 p ^ Q ^ 9 i+ V 2 ,j- V2+P.k- V2+q ~ 6 i" \ ,j- V2+P,k- V^+q- 


^ G x 9) i,j,k 


(A. 7) 


1 1 


" Ah q | Q p ^ 0 ^ 9 i- V2+q,j+ V 2 »k- V 2 +P 9 i- V 2 +q,j- l h »k- V 2 +P^ 


- (G h 0) . . 


(A. 8) 


1 1 


3 1 * * 

i,j ,k " Ah p ^ Q ^ 9 i- V2+P,j“ V2+q>k+ V 2 9 i- V2+P»j“ V2+q,k- V 2 ^ 


= ^ 9) i,j,k* 


(A. 9) 


Then (A. 4) implies that 0 satisfies the following discrete Poisson equation 


A. 0 - (D G h 0 + D G h 0 + D G h 0) . = <{>. . 

h xx y y z z i,j,k i,j,k 


where A^ is the discrete Laplacian shown in Figure 2. (a^ is the same as 

the "one-point quadrature stencil" described in [15].) 

The equation above only involves 0 at the centers of boxes. Since all 

boundary conditions are natural, we can extend the grid to cover £2 and let 

0, . i, . i, , be zero for points outside 0. We now show that such a 

i+ V2 » j+ V2 »k+ V2 

0 exists. 


A^ t|> = f in £2 and \p ■ 0 outside £2. 


(A. 10) 


[Here f is constant over each box Then, multiplying both side by 

t|> and performing a discrete integration by parts, we obtain 


- 21 - 


-{ l 

(i,j,k)en 


( G h 9)J . + 

x 'i,j,k 


(c h 0) 2 . + 

y i»j»k 


(G * ■ l (£6, i+ V 2 ,J+ V 2 ‘/2 


Put f = 0. Then, since 0 = 0 outside 0, by starting from a corner of fl 
and working through 0, we see that 0=0 in 0, Hence, the Fredholm 

Alternative, (A. 10) has a unique solution for each f. We therefore obtain a 
solution for (A. 4) . Defining v*, v^, v^ by (A.7)-(A.9) we get the standard 
estimate 

3 

h 2 l l | 2 <CH<f>ll 2 (A. 11) 

£=1 i,j,k 


for some absolute constant 0 < C < °°. In addition, defining v at the 
centers of boxes by (A. 6), we obtain 


3 

l 

£=1 


A 


(i+ V2 , j+ V2 »k+ V2 ) e a 


I 1 

|v i+ 1/2 ,j+ V2 >k+ V2 


| 2 < c 2 {n<f>ii 2 1 


+ 



Letting v be the continuous piecewise linear function whose value at each 
node is defined as above, we get 


3 

l 


1=1 


ivh 2 < 



and 


div v 




Finally, using the inverse inequality for the uniform grid on 0, we obtain 
for <{> £ S h , 

n<H 0 < ch - 1 n<{.ii-i. 


This proves the theorem 
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Figure A.l 
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Figure A. 2 
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